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The inspiral of a stellar compact object into a massive black hole, an extreme-mass-ratio inspiral, 
is one of the main sources of gravitational waves for the future space-based Laser Interferometer 
Space Antenna. We expect to be able to detect and analyze many cycles of these slowly inspiraling 
systems, which makes them truly high precision tools for gravitational-wave astronomy. To that 
end, the use of very precise theoretical waveform templates in the data analysis is required. To build 
them we need to have a deep understanding of the gravitational backreaction mechanism responsible 
for the inspiral. The self-force approach describes the inspiral as the action of a local force that can 
be obtained from the regularization of the perturbations created by the stellar compact object on 
the massive black hole geometry. In this paper we extend a new time-domain technique for the com- 
putation of the self-force from the circular case to the case of eccentric orbits around a non-rotating 
black hole. The main idea behind our scheme is to use a multidomain framework in which the 
small compact object, described as a particle, is located at the interface between two subdomains. 
Then, the equations at each subdomain are homogeneous wave-type equations, without distribu- 
tional sources. In this particle-without-particle formulation, the solution of the equations is smooth 
enough to provide good convergence properties for the numerical computations. This formulation is 
implemented by using a pseudospectral collocation method for the spatial discretization, combined 
with a Runge Kutta algorithm for the time evolution. We present results from several simulations of 
eccentric orbits in the case of a scalar charged particle around a Schwarzschild black hole, an excel- 
lent testbed model for testing the techniques for self-force computations. In particular, we show the 
convergence of the method and its ability to resolve the field and its derivatives across the particle 
location. Finally, we provide numerical values of the self-force for different orbital parameters. 
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I. INTRODUCTION 

One of the main subjects of activity in the emergent 
area of Gravitational Wave Astronomy is the physical 
description of the different astrophysical and cosmologi- 
cal sources of gravitational radiation that are expected 
to be detected by the present generation (LIGO [T], 
GEO600 0, VIRGO 0) and future/planned (LCGT g], 
LISA ET [5], etc.) laser interferometer observato- 
ries. The description of sources of gravitational waves 
plays an important role for several reasons. On the one 
hand, we need to understand well these sources in or- 
der to convert the information from their observations 
in gravitational waves (and perhaps in some electromag- 
netic bands) into new scientific knowledge. On the other 
hand, in order to detect and obtain accurate physical in- 
formation from these observations we need to use data 
analysis techniques, like the matched filter, that require 
the use of precise waveform templates in order to separate 
the signals from the noise and to estimate appropriately 
the source parameters. 

In the case of the future Laser Interferometer Space 
Antenna (LISA) [5] (see [7J |S] for the proceedings of 
the 7th international LISA symposium), the main ex- 
pected sources of gravitational waves are: Galactic com- 



pact binaries, mergers of Massive Black Holes (MBHs), 
the inspirals of stellar-mass compact objects into MBHs, 
and stochastic gravitational-wave backgrounds of diverse 
origin. This paper is devoted to the description of the 
dynamics of the third type of source which, due to 
the extreme mass ratios involved, are usually known as 
Extreme-Mass-Ratio Inspirals (EMRIs) . The EMRIs of 
interest for LISA are associated with relaxed galactic nu- 
clei that have a dense stellar population (a stellar cusp 
or core) around a hypothetical central MBH. In this sit- 
uation, there is a non-negligible probability for a stellar- 
mass compact object (SCO), with masses in the range 
m ^ 1 — 50Af Q , to end up in a slow inspiral towards 
the MBH (see also (9] for details and references on the 
astrophysics of EMRIs). The range of the MBH mass 
of interest for LISA is M ~ lfj 4 — 10 7 M Q , which means 
that the mass ratios involved in EMRIs are in the interval 
fx = m/M ~ 1CT 7 - 1CT 3 . 

In brief, the situation we face in trying to detect EMRIs 
with LISA is the following: They are systems that dur- 
ing the last year of inspiral before plunge complete of the 
order of 10 5 cycles [10J, so they fall into the ~ 0.1 mHz 
part of the LISA band. Therefore, the mass ratio im- 
plies not only the existence of very different spatial scales 
(e.g. the size of the SCO as compared to the size of the 
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MBH) , but also of two very different time scales, namely 
the orbital time scale versus the inspiral time scale, so 
that T orbitgl /T. irgi ~ \i = m/M. Despite the difficulties 
introduced by the extreme mass ratios involved in the 
problem, there is an advantage to it: we can assume that 
the spacetime geometry is the one due to the MBH with 
deviations produced by the SCO. Then, one can use per- 
turbation theory to deal with EMRIs. In this context, the 
inspiral can be described as the action of a local force, the 
so-called self-force, according to the equation of motion 
derived in [111 IT2"] . known as the MiSaTaQuWa equation 
of motion. The MiSaTaQuWa equation can be used as 
the foundations of a self-consistent scheme to simulate 
an EMRI event in an adiabatic way by coupling it to the 
Partial Differential Equations (PDEs) describing the per- 
turbations produced by the SCO (together with a scheme 
to identify the regular part of the perturbations respon- 
sible for the inspiral). For recent studies on these issues 
see [13 15 . For reviews on different aspects of EMRIs 
see [MIB]- 

Given the amount of cycles required for EMRIs and 
their computational cost we cannot expect to generate 
complete waveform template banks by means of full self- 
force calculations. Instead, the goal of these studies 
should be to understand all the details of the structure of 
the self-force so that we can formulate efficient and pre- 
cise algorithms to create the waveforms needed for LISA 
data analysis, perhaps complementing present approxi- 
mate schemes [TT?1 - |2"2"] . 

Moreover, due to the large number of cycles that will be 
available in future LISA EMRI observations, we expect 
to determine the EMRI physical parameters with high 
precision [TU]. Given also that the expected event rate is 
in the range 10 - 10 3 EMRI/yr (HH |H] up to a redshift 
z w 1, these observation will have an important impact 
not only on astrophysics (see [HI and references therein), 
but also on cosmology and fundamental physics [25 , from 
precise measurements of cosmological parameters |261 |2"T] 
to tests of the geometry of MBHs [28H34] and/or the 
gravitational interaction [351 I36| . 

If Intermediate-Mass Black Holes (IMBHs) (37] d] ex- 
ist (with masses in the range 10 2 — 1O 4 M ), as sev- 
eral recent observations suggest, they can form binaries 
with intermediate mass ratios resulting in the so-called 
Intermediate-Mass-Ratio Inspirals (IMRIs). LISA will 
be sensitive to one type of IMRIs: That an IMBH into a 
MBH. In contrast, the inspiral of an SCO into an IMBH, 
another type of IMRIs outside the LISA band, is of in- 
terest for the next generations of ground-based detec- 
tors [391 3D]. We can obtain an approximate descrip- 
tion (probably not accurate enough for data analysis pur- 
poses) of IMRIs by using the same methods that are used 
for EMRIs. 

In the calculation of the self-force, it is reasonable to 
neglect the structure of the SCO and describe it as a 
point-like object. This means that the retarded metric 
pertubations generated by the SCO diverge at the parti- 
cle location, and hence they need to be regularized. An 



algorithm that has been widely employed is the mode 
sum regularization scheme |41ll4"8"] . It substracts, mode 
by mode (harmonic £-modes) from the retarded field, the 
singular part of the perturbations that does not con- 
tribute to the self-force. We therefore need a method for 
computing the different multipoles of the retarded part 
of the metric perturbations induced by the SCO. Due to 
the complexity of these linear equations, solutions can 
only be obtained numerically, either using frequency- or 
time-domain methods. Both have advantages and dis- 
advantages. Because EMRIs are expected to have high 
eccentricities and frequency-domain methods have shown 
difficulties in this case, we concentrate on time-domain 
methods. Time-domain methods have difficulties related 
to the resolution of the SCO (see, e.g. [IS]), since many 
of them need to introduce an artificial scale in order to 
model the distributional energy-momentum tensor that 
describes the SCO. 

To solve the problems related to the artificial scale as- 
sociated with the SCO some new techniques have been 
recently proposed [50- 53] which are similar to the punc- 
ture method introduced in numerical relativity |54H56| . 
However, they do not completely solve the problem of the 
spatial resolution associated with the source terms. We 
recently introduced |57[ [58] a new technique, for the case 
of circular orbits, that completely avoids the introduction 
of a spatial scale. The method, described in detail in [55] 
(henceforth Paper I), is based on the splitting of the spa- 
tial computational domain (the one-dimensional domain 
associated with the radial coordinate) into several sub- 
domains, in such a way that the SCO is located at the 
interface of two such subdomains. For the circular case, 
the SCO is at rest in the spatial domain. The outcome 
of this setup is that the equations at each subdomain are 
homogeneous wave-like equations. In this way, we ensure 
the smoothness of the solutions in each domain and at the 
same time we avoid artificial scales associated with the 
SCO. The SCO appears in the boundary/matching con- 
ditions between the subdomains that surround it. This 
setup is ideal to use high-precision numerical methods, 
which normally require a high degree of smoothness in 
the solution. In Paper I we implemented these ideas 
for the case of a charged scalar particle in circular or- 
bits around a non-rotating MBH using the PseudoSpec- 
tral Collocation (PSC) method for the spatial discretiza- 
tion (a similar technique has also been proposed recently 
in [53]). In this sense, the smoothness of the solutions of 
the homogeneous wave-like equations involved, ensures 
the exponential convergence of the method, as shown in 
Paper I. This study has further demonstrated the high 
precision of the method by comparing results with previ- 
ous calculations of the self-force on a charged particle in 
circular motion [481 RH"] . Finally, it was shown that 
by adjusting the subdomains and the number of collo- 
cation points, one also obtains an efficient method (in 
terms of computational time) for the calculation of the 
self-force. These techniques can be directly imported to 
the gravitational (and electromagnetic) case. 
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The main aim of the present paper is to show that 
the techniques shown in Paper I can be extended to the 
case of generic eccentric orbits. The key ingredient of our 
formulation, which we call the Particle-without-Particle 
formulation, is to divide the computational domain into 
subdomains and to locate the particle at the interface of 
two of them in such a way that the equations do not have 
a distributional term that deteriorates the smoothness of 
the solution. The challenge is to ensure communication 
between the domains in an appropriate way to take into 
account the particle effects in a precise way. The main 
challenge in comparison with Paper I is to adapt these 
techniques to the case where the particle is no longer 
at rest in the spatial domain, as its radial coordinate is 
changing with time. This will be achieved by adapting 
the mapping between the physical domain (correspond- 
ing to the radial coordinate) and the spectral domain 
(corresponding to the coordinate that labels the collo- 
cation points used in the PSC method) such that the 
particle is moving in the physical domain but it is at rest 
in the spectral domain. In this way, we can apply all the 
machinery shown in Paper I. Here we present the ingre- 
dients of the new formulation in order to handle eccentric 
orbits: the equations, the junction conditions, different 
ways of implementing it numerically, etc. In particular 
we show that the spectral convergence is preserved. We 
also show that we can evolve the system in time with 
long-term stability and a good resolution of the jumps in 
the derivatives of the solution induced by the presence of 
the particle. In particular, we provide values of the self- 
force components, at the pericenter radius, for different 
eccentric orbits. 



The organization of this paper is as follows: In Sec. |TT] 
we introduce the formulation of the physical model, a 
charged scalar particle orbiting a Schwarzschild MBH, 
including the basic formulae for the computation of the 
self-force via the mode-sum regularization scheme. In 
Sec. IIIII we introduce the basic mathematical formula- 
tion of our method, based on the division of the physical 
domain into subdomains and on the communication be- 
tween them. In Sec. |IV| we describe how to implement 
the mathematical formulation of Sec. |III| using the PSC 
method. In Sec. [V]we show the performance of the nu- 
merical codes that we have designed to implement the 
new scheme as well as results from the computation of 
the self- force. We conclude with a summary and discus- 
sion in Sec. |VI[ including future work. Three appendices 
complete this paper: In Appendix |A"| we summarize the 
main equations that describe the geodesic motion that 
we use in this work. In Appendix [B] we provide the main 
formulae related to the structure of the singular field 
and the associated regularization parameters. Finally, 
in Appendix [C] we just provide the expressions for the 
time derivatives of the jumps in the characteristic vari- 
ables. Throughout this paper we use the metric signature 
( — ,+,+,+) and geometric units in which G = c = 1. 



II. CHARGED SCALAR PARTICLE ORBITING 
A SCHWARZSCHILD MBH 

In order to simplify the presentation of our techniques 
we use a simplified EMM model, corresponding to a 
charged scalar particle orbiting a non-rotating black hole. 
In this model the SCO is represented as a particle, with 
charge q associated to a scalar field $, and the MBH 
is described by a fixed spacetime geometry, in the sense 
that it is not affected by the charged particle. Since we 
are interested in the dynamics of EMRIs in General Rel- 
ativity, the geometry of the non-rotating MBH is given 
by the Schwarzschild spacetime: 

ds 2 = f(-dt 2 + dr* 2 )+r 2 dn\ f = l- — , (1) 

r 

where M is the MBH mass, dfl 2 = dQ 2 + sin 2 9 dip 2 de- 
notes the metric of the two-sphere, and 

r* =r + 2M\n(— - l) , (2) 

is the so-called radial tortoise coordinate. The particle 
induces the scalar field <E> satisfying the following wave- 
like equation (see, e.g. |17|): 

g a ? V Q V^$(x) = -47T9 f dr 6 4 (x, z(t)) , (3) 

where g°^ denotes the Schwarzschild inverse metric |l]) 
and the associated canonical connection; r denotes 
proper time associated with the particle along its time- 
like worldline 7, with coordinates x 11 — z^(t) and 
8±(x, x r ) is the invariant Dirac biscalar distributional 
functional in Schwarzschild spacetime, defined as 

(<5 4 , /) (x) ee J d 4 x'y/^)f(x')S 4 (x', x) = f{x) , (4) 

with g being the metric determinant. 

The field equation Q has to be complemented with 
the equation of motion for the scalar charged particle: 

dv^ dz^ 
m^T =F» = gGr+u"u") (V„*)| 7 , «* = ~, (5) 

where m and are the mass and 4-velocity of the par- 
ticle, respectively. The coupled set of equations formed 
by the Partial Differential Equation (PDE) for the scalar 
field (J3J and the Ordinary Differential Equation (ODE) 
for the particle trajectory (|5| constitute our testbed 
model for an EMRI. The SCO generates scalar field ac- 
cording to pj, which in turn, affects to the SCO motion 
according to (pj, that is, through the action of a local 
force. This mechanism is the analogous of the gravita- 
tional backreaction mechanism that produces the inspiral 
via the gravitational self- force. 

However, this is not the end of the story since the re- 
tarded solution of ([3]) is singular at the particle location, 
whereas the force in Eq. (I5J) involves the gradient of the 
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field evaluated at the particle location. Therefore, as they 
stand, Eqs. ^ and (5]) are formal equations that require 
an appropriate regularization to become fully meaning- 
ful. Following [62 , the retarded field can be split into two 
parts: a singular piece, $ s , which contains the singular 
structure of the field and satisfies the same wave equation 
as the retarded field, i.e. Eq. (|3]), and a regular part, $ R , 
that satisfies the homogeneous equation associated with 
Eq. pj. As it turns out, <1> R is regular and differentiable 
at the particle's position and is the solely responsible of 
the scalar self- force 162 . We can therefore write 



(0) 



which gives a definite sense to the equations of motion of 
the system. 



In solving the field equation 
the spherical symmetry of the MB 
the solution in scalar spherical harmonics (see Paper I 
for conventions): 



we can make use of 
1 geometry to expand 



K = EE & m (t,r)Y lm {6, V ) . 



(7) 



1=0 m=-e 



The advantage of this expansion is that each harmonic 
mode, $ (t, r), is decoupled from the rest and satisfies 
the following 1 + 1 wave-type equation: 



If 



d 2 



dt 2 dr* 
where 



+ 7^2 " ) = s s ( r ~ sW) . ( g ) 



V t {r) = f(r) 



£(£ + 1) + 2M 



(9) 



is the Regge- Wheeler potential for scalar fields on the 
Schwarzschild geometry and 



aim _ 4n( lf(r p ) vtm( TT 

6 " ( 2'^ J 



(10) 



is the coefficient of the singular source term |76| due to 
the presence of the particle, and the bar denotes com- 
plex conjugation. Here we have assumed, without loss 
of generality, that the particle's orbit takes place in the 
equatorial plane 9 = tt/2. Moreover, r p and tp p denote 
the radial and azimuthal coordinates of the particle, and 
are functions of the coordinate time t. 

The expansion in spherical harmonics is useful not only 
for simplifying the construction of the retarded solution, 
but also in terms of constructing its regular part, $ R . 
Indeed, it turns out that each harmonic mode of the re- 
tarded field, <l> ferl (t,r), is finite and continuous at the 
location of the particle. It is the sum over £ what di- 
verges there. Here is where the mode-sum regularization 
scheme |41[ |45j [46] comes into play, as it provides ana- 
lytic expressions for the singular part of each £-mode of 
the retarded field. These expressions for the singular field 



are valid only near the particle location. The regularized 
self-force is thus obtained by computing numerically each 
harmonic mode of the self-force and subtracting the sin- 
gular part provided by the mode-sum scheme. Thus, the 
regular part of the gradient of the field, V a $ R = $ R , is 
given by 

oo 

$ r (^(t)) = lim IKM " a&V)) ■ (11) 
Introducing, 

£ 

= E V a (* tm (t,r)Y tm (6,<p)), (12) 
m=-i 

the structure of the singular field can can be written as: 



+ 



£+- I .1.. -B. 



£+\ 



(£-m+D 



(13) 



The expressions for the regularization parameters, A a , 
B a , C a , and D a , can be found in the literature for generic 
orbits HHJ E3J EJ. They do not depend on £ but de- 
pend on the trajectory of the particle. The three first 
coefficients of ( 13 1 are responsible for the divergences, 



whereas the remaining terms converge to zero once they 
are summed over £. The expressions for the regulariza- 
tion parameters that we use in this paper are listed in 
Appendix [B] 



III. THE PARTICLE- WITHOUT-PARTICLE 
FORMULATION 

Before we introduce the multi-domain framework for 
the numerical calculations that we present in this paper, 
it is important first to look at its mathematical founda- 
tions. To that end, we only need to consider two subdo- 
mains of the physical domain, that is 



r* e (— oo, oo) = (—oo, r*) U (r*, +oo) 



(14) 



where r * — r* (t) is the radial tortoise coordinate of the 
particle. This particular division obeys to the fact (see 
Paper I) that the key point of our technique consists in 
locating the particle between two subdomains, which au- 
tomatically makes all the wave-type equations governing 
the dynamics of the scalar field modes to be homoge- 
neous. As a direct consequence, we avoid all the prob- 
lems associated with source terms with support at just 
one point and prevents us from introducing an artificial 
scale associated with the charged particle. 

The problem is thus reduced to the junction condi- 
tions for the field across the particle. It is well-known 
(see, e.g. |65| ) that discontinuities in hyperbolic equations 
(like the wave equation) can only appear and propagate 
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along characteristics. In order to analyze the junction 
conditions, it is therefore convenient to adopt a hyper- 
bolic formulation of the equations for the scalar field. To 
that end, we introduce the variables 

U=(fl> lm ,<tP n ,<fP n ) = (r¥ m ,d t ^ em ,d r ,ip tm ). (15) 

at each of the two subdomains in ( fl4| , these satisfy the 
first-order system of PDEs (see [66 for an analysis in 
second-order form): 



d t U = A-d r M + M-U + S, 
where the matrices A and B are given by 



(16) 














1 






-G 





■)■ 


B = 









(17) 




l 






w 










and the vector S is 



S= 0,- 



s ( 



f(r p 



-6(r*-r*Jt)),0 



(18) 



One can check that the system of PDEs given in Eq. ( 16 ) 
constitutes a first-order symmetric hyperbolic system. In 
order to analyze the jumps across the particle in terms 



of the variables of Eq. (15 1, we introduce the following 



splitting based on the domain division of Eq. ( 14) 

U(ty) = U_(t,r*)G(r* p (t)-r*) 
+ U + (t,r*)Q(r*-r* p (t)), 



(19) 



where U_ and U , are going to be smooth functions on 
(— oo,r*) and on (r*,+oo) respectively, and O denotes 
the Heaviside step function. Introducing this ansazt 

and using the following 



[Eq. (|TQ[)] into the system (16 



definition of the jump in an ar 



□itrary quantity A 



[X] p = li m< A + (*,r*)- lim,A_(t,r*) 

r — >r p r -^r p 

we obtain the evolution equations for U±: 
d t U ± = A-d r ,U ± +B-U 



± ■ 



(20) 



(21) 



and the expressions for the jumps in U across the particle 
(see also [66 and Paper I): 



W m ] P = o, 

L h (l 



{V 1 ™] = 77 



(22) 
(23) 

(24) 



where the dot denotes differentiation with respect to 
the time coordinate t, /„ = f(r p ) and S is given in 
Eq. (10 1. As we can see, the field is continuous across 



the field have a jump, which depend on the particle's 
speed and satisfy the following advection-like equation: 



[< 



■] =0. 



(25) 



For r* — const., i.e. r* 



0, we recover the jump condi- 
tions for the circular case studied in Paper I. 



As mentioned above, the jump conditions (22l-(24l 
have to be imposed on the characteristic fields because 
discontinues propagate along the characteristics. In our 
case, the characteristic fields are: i// m (whose charac- 
teristic surfaces are the spacelike surfaces {t = const.}), 



U = <fi tm — (p lm (whose characteristic surfaces are the 
null surfaces {t — r* = const.}), and V lm — (j) im + ip tm 
(whose characteristic surfaces are the null surfaces {t + 
r* = const.}). See an illustration of this in Figure [2] 

An alternative description of the system, which we 
shall also use in this work, is to choose a set of vari- 
ables associated with the characteristic fields. The obvi- 
ous choice is 



N = (t/j tm ,U im ,V* m ) 



(26) 



As before, we can apply the same splitting of Eq. ( 19 1 to 
these variables and follow the same procedure to obtain 
a set of equations for N±: 



d t N ± = C-d r ,N ± 
where the matrices C and D are 



■ N 



± i 














-1 


■ 


D= ( 












1/2 1/2^ 
-V e ( 
-V, ( 



and the junction condition of Eq. ( 22 ) and 



\u im ] = 

\v lm ] = p 



(27) 
(28) 

(29) 
(30) 



We can now establish what we call the Particle with- 
out Particle (PwP) formulation of the problem. By split- 
ting the physical domain from the particle location as in 
Eq. (14) we can introduce in a natural way the split- 



ting in the dynamical variables given by Eq. ( 19 1 (and 
the equivalent one for the characteristic variables TV). In 
this way, the restriction of the global variables to the 
subdomains to the left of the particle, (— oo,r*), and 



to the right of the particle, (r*,+oo), U + {N + ) and 
U (N_) respectively, satisfy homogeneous hyperbolic 
equations (21 1 [Eq. (27 1 for the characteristic variables]. 



This means that we have got rid of the distributional 
source terms, due to the particle, that appeared in the 
global equations (16 1. Then, the presence of the particle 



is introduced in the communication between the variables 
of different subdomains through the junction conditions 
(|22]|-p4| [Eqs. ((22l, p9l), and Q for the characteristic 



the particle, and only the time and radial derivatives of variables]. 
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IV. A PSC METHOD IMPLEMENTATION OF 
THE PWP SCHEME 

We now introduce a multidomain computational 
method that implements the PwP scheme described in 
the previous section. To that end, we follow the dis- 
cussion of Paper I. generalizing the methods introduced 
there to include the case of generic orbits, in which the 
particle is not at rest in the physical domain. 

The first step is the truncation of the physical domain 
(as opposed to compactification) . Instead of consider- 
ing the domain ( 14 1 we consider the following domain: 
il = [r*,r*], where r* corresponds to the truncation in 
the direction to the MBH horizon (at r* —> — oo) and 
r* corresponds to the truncation in the direction to spa- 
tial infinity (r* — > +00). Next, we divide this truncated 
domain into a number D of subdomains: 



D 



n 



A,L J ' A : R| 



(31) 



where r * L and r* R are the left and right boundaries of 
the subdomain f2 A . Following the PwP scheme we place 
the particle at the interface between two of the subdo- 
mains. Now, let us assume that the particle is located 
at the interface between the subdomains f2 p and f2 P+1 , 
then 

(32) 

The main difference with respect to the circular case of 
Paper I is that now r* is a function of the time t. This im- 

must be functions of t as well, 



P+1,L 



plies that r* R and 
and the subdomains will in general change with time. 
We do not need that all subdomains change with time 
however, so we distinguish between dynamical and non- 
dynamical subdomains according to whether or not any 
of their boundary nodes moves with time. This multido- 
main framework is illustrated in Fig.[l] In the remainder 
of this work we restrict ourselves to the case in which 
only the two subdomains that surround the particle are 
dynamical, i.e. only r* R and r* +1 L are time dependent. 

The first ingredient in applying the PSC method is 
the discretization of the spatial domain (see [67J for a 
multidomain PSC method for elliptic equations). In our 
scheme, each subdomain is discretized independently, by 
using a Lobatto-Chebyshev grid (see Paper I for details on 
our particular implementation, and e.g. [68 for general 
details on the method) , with collocation points given by 



7T I 



(i = 0,...,N). 



(33) 



In order to complete the discretization of the physical do- 
main, we need to specify how the spectral coordinate X, 
associated with the Chebyshev polynomials, is related 
to the physical coordinate r* . This is the key ingredi- 
ent in our scheme; it is here where we establish how to 
keep the particle at the interface of two subdomains even 
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FIG. 1: The figure shows the structure of the one-dimensional 
spatial grid, the division in subdomains and the location of 
the particle at the interface between two of them, rp(t). In 
a generic orbit, the particle moves in a periodic way between 
the minimum (pericenter, denoted by r peri ) and the maxi- 
mum (apocenter, denoted by r apo ). In this example, only the 
domains on the left and on the right of the particle are dy- 
namical, i.e. they are the only ones whose coordinate size 
changes in time. 



when it is moving in the radial direction. This informa- 
tion is encoded in the mapping from the spectral and 
physical domains. We have effectively traded particle 
motion for coordinate mappings. Given that the particle 
will be in general moving, it is convenient to prescribe a 
spacetime mapping. Since we discretize each such sub- 
domain independently, and the spectral domains are al- 
ways the interval [—1, 1], there are going to be as many 
mappings as subdomains. Then, for a given subdomain 
f2 A = [r* L , r* R l (A — 0, . . . , D), the spacetime mapping 
to the spectral domain is taken to be linear and given by: 



X A : Tx [rJ iL ,rI iB ] 



(t,r* 



T x [-1,1] 
(T,X A ) 



(34) 



where T = [t ,tj] is the evolution time interval and 

T(t) = t , 
X A (t,r*) 



2r* - r A , L - r AiR 



(35) 



Here, the time dependence in X A comes from the varia- 
tion of r* L or r* R in time. The inverse (linear) mappings 
from the spectral domain to each of the subdomains fl A 
are given by 



r*L : Tx [-1,1] ^Tx [r AiL ,r AiR ] 



(T,X) 



{t,r* 



(36) 
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and 



t(T) = T , 



A.R ' A.L 



X 



(37) 



X A ■ is the i-th collocation point at the domain f2 A , given 



by Eq. (33 I, and 



(d r ,X A 



(45) 



In the PSC method we have two expansions (represen- 
tations) of our variables: The spectral and the physical 
representations. The former expands the variables in a 
basis of Chebyshev polynomials, {T n (X)} (X G [—1,1], 
n = 0, . . . ,N). Thus, our variables U can be approxi- 
mated by 



A' 



U N (t,r*) = J2"n(t)T n (X A (t,r*)), 



(38) 



where the a n are the spectral coefficients which, in our 
case, will be time dependent. The physical representation 
approximates the solution as 



U N (t,r*) 



N 

£ 

i=0 



(39) 



where C^X) are the cardinal functions associated with 
our choice of collocation grid (Lobatto-Chebyshev) which 
satisfy the relation C^Xj) = 5^. In this representation, 
the coefficients {t^} i=0 N are the values of our vari- 
ables 

Uit^iTitlX^^US) (i = 0,...,N), (40) 

at the collocation points. 

In these two representations we have two types of un- 
knowns, the coefficients {a n } n=a N and the collocation 
values {L/j} i=0 N . These N + 1 unknowns are fully 
determined, at each subdomain, by imposing our set of 
PDEs (21 I [or Eq. (27)] at each of the collocation points. 

Then, the discretized version of the evolution equations 



(21 1 at each subdomain f2 A takes the following form: 



(d x U)i 



u, 



(41) 



where there is no summation over i and the matrices 
and Qj are given by 



/0 

o -(d t x A \ [d r ,x A \ 
\0 (d r ,X A ) i -(d t x A ) iy 



-V^ 1 
V 



(42) 



(43) 



the quantity V^' 1 is given by 

V i A > i (t) = Vz(r*(T(t),X A M, 



(44) 



S A ' T 



(46) 



\-\ j- j 1 il r p 



A,L 7 

otherwise, 



7\ 



1 if r* - . 

'p ' A,R ' 

otherwise. 



For a given set of initial data 

U° = U(t = t ,r*(T,X A J) 



(47) 



(48) 



the system of N + l ODEs defined by Eq. (41 I determines 



uniquely the values of our variables U at the collocation 
points, i.e. the /V+l unknowns, which in turn determines 
uniquely our approximation to the solution via Eq. ( |39[ ) . 
Looking at the expressions for the matrices P i and Q i 
[Eqs. (42) and (43 1] we see that their components are ex- 
plicitly time dependent in the case of a dynamical subdo- 
main (in the case of this paper, only those that surround 
the particle). 

A similar system of equations can be obtained for the 
characteristic variables TV [Eq. (26l], 



(49) 



where again there is no summation over i and the matri- 
ces KL and S„- are 





(d t X A \ - {d r *X A \ 







■{d t X K \ + (d r *X A ). 



-A,i 








(3t* A )i 

{s r »x A ). 







(51) 



/ 



However, the ODEs given in Eq. (41) do not const! 



tute a complete description of our system. We still need 
to incorporate all the necessary boundary conditions to 
adequately describe the system. First of all, we have 
the global boundary conditions, at r* and r* to pre- 
vent incoming signals from outside the physical domain, 
n = [?'* , r*] , sometimes known as absorbing boundary 
conditions. In this work we try to push the boundaries 
as far as possible (typically in such a way that, given the 
evolution time, tf — t Q1 the boundaries are not in causal 
contact with the particle) and impose Sommerfeld out- 
going boundary conditions 



0, 
0. 



(52) 
(53) 
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These are exact boundary conditions only when V e = 
or r* — > —oo and r* — > +oo. 

The second type of boundary conditions that need to 
be imposed are the junction conditions between the dif- 
ferent subdomains. In this paper we present two differ- 
ent methods of imposing them: (i) The penalty method, 
which has already been used in Paper I, and (ii) direct 
communication of the characteristic fields. In the remain- 
der of this section we describe both methods and the 
modifications that have to be included for their imple- 
mentation in the system of ODEs (41 1. We recall that 



in this paper we consider the case in which the only dy- 
namical subdomains are those two surrounding the par- 
ticle, fl p and f2 P+1 . In this way the particle is located at 
their interface as described by Eq. (32 1. We also have: 



r* N and 



' a+i,o- 



A. The Penalty method 

The basic idea behind the penalty method is to mod- 
ify the equations by adding terms that dynamically drive 
the system to a state in which the junction conditions are 
satisfied to a high degree of precision. These additional 
terms are proportional to the conditions that one is in- 
teresting in imposing. In our case they are proportional 
to the junction conditions. Details of the method can be 
found in Paper I (for generalities see, e.g. [69 ). 

An important point in using the penalty method is 
that since we are dealing with hyperbolic equations, the 
junction conditions have to be imposed across the char- 
acteristics (see Figure H| , and this determines which com- 
binations of the junction conditions have to be used in 
each equation. 



t + r* 



k 



const. 

\ v = 4> + <p 



Time (t) 




t — r = const. ,-' 

,<\U / 

/\U / 
/ U = cj)-<p 



Particle Nodes 



r, 



FIG. 2: Characteristic structure of our field equations. The 
picture shows the characteristics (t ± r* = const.) of the hy- 
perbolic equations that describe the system [Eq. \16\ ], the 
characteristic fields (i//™ , U lm = <f) lm - </ m , and V tm = 
4> lm + <p em ), and their propagation direction. 

In what follows we list the different equations for the 
different cases that appear depending on whether we are 



dealing with a subdomain to the left or to the right of 
the particle, and also on whether the variables are associ- 
ated with boundary points or not. That is, we have four 
basic cases (other possibilities can be derived from these 
cases): (i) Inner points for the subdomain to the left of 
the particle, i.e. i7 P ; (ii) inner points for a subdomain 
to the right of the particle, i.e. f2 P +i; (iii) left boundary 
point; and (iv) right boundary point. 

(i) The equations for inner points (i = 1, . . . , N — 1) 
valid at the dynamical subdomain at the left of the par- 
ticle (to simplify the notation we have dropped the har- 
monic indices t and to) are: 



r 

p 2 

r 

p 2 

1 + A D 



(dr*¥>p)i + (0r-<£p)i 



(54) 

(55) 
(56) 



where for any quantity A we compute the spatial deriva- 
tives as 



(d r .A P ) i = {d r .x J .) i (d x A P ) l , 



(57) 



and (d r ,X p ) i is given in Eq. (45) 



(ii) Equations for inner points at the dynamical sub- 
domain at the right of the particle: 



d T ip p+li 
d T (t> P 



l-X 



r 



p+i, 



2 

1 — Xp+ij 
2 

A-X v . 



Vp+m + 0p+i,» > ( 58 ) 
(9 r ^ p+1 ). + (9 r ,^ p+1 ). 

(59) 

(d r »(p p+1 ) i + (d r *<t> p+1 ){6Q) 



The equations at the inner points of any non-dynamical 
subdomain are obtained by setting r* = either at the 
set of equations ( 54 1-( 56 1 or at the set of equations (58)- 



(iii) Equations valid at the right boundary of the dy- 
namical subdomain at the left of the particle, that is, at 
r* =r P , N = r* p (t): 



<9 T V F 



T 4,' R (^p,n - ^p+i,o) > 
r* (9 r ,0 P ) N + (d r ,ip P ) N 



(61) 



P,R 

p,p+i 



[<t>] 



lP.P+l 



+ 1,0 + ^p + i.o) 



V 



i<p. 



(62) 



<9 T ¥>P,N = 



r* (d r ,ip p ) N + (d r ,(j) p ) N 

»,N + ^P : N - (0P + 1.O + <Pp+l,o) 



2 



p,p+i 



(63) 
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,12 , and t^ ,r , are penalty (constant) coeffi- 
cients associated with the evolution of tp p N (T), 4> p N (T), 
and <p PN (T) respectively. The quantities [0] P ,' P+1 and 
[<p]p ,P+1 are the jumps across the subdomains f2 P and 



flp+i, given in Eqs. (23 I and (24 1 



(iv) Equations valid at the left boundary of the dy- 
namical subdomain at the right of the particle, i.e. at 

r — r p+i,o — r p\ l y- 



d T ip P+10 — r p fp+ifi + P +i,o 



^p+i,o 



P+1,L / I 

r* (d r »<f> 



p+i,o 



p+Uo n 

Wl,0 

+ [<p] 



■ ^p,n) , 

- - (4 



T/P + 1,0 , 

v, V 



d T ip 



p+1,0 



[4>}p . L^jp 

r* (d r »tp P+ i) + (d r ,(j} p+1 ) 

T P+1,L 

— [0p+i.o - Pp+i,o - G 



p,p+l 



1P,P+1 
p 



(64) 

p+1,0 

(65) 
(66) 

(67) 



where again t1 + ' l , t ( p+1 ' l . and t£ +1,l are penalty co- 
efficients associated with the evolution of ip P+1 o(^)> 
P+1 (T), and v? P+1 oC? 1 ) respectively. At this point we 
comment on the equation for ip im . Here, we are enforc- 
ing the continuity of this variable via a penalty term. 
However, as we already mentioned in Paper I, there is an 
alternative to this, which consists in taking a sort of av- 
erage of the right-hand side of the equation for ip im over 
the two subdomains that we are trying to communicate 
or match. The alternative equation for ip tm would look 
as follows: 



d T tp 



1 



p + 1,0 



~jp (</Vn + 9 



p+1,0, 



+ 



1 



(68) 

and the same applies to the equation for ip P . In this way, 
-0 is continuous, to machine precision, by construction. 

Finally, the outgoing (global) boundary conditions de- 
scribed by Eqs. (52) and (53 I are equivalent to the fol- 
lowing equations: At the horizon r* = r* : 



d T 4>i, H = d T ip liH , 
9 T iPi, H = (<9 r *<£i) H > 

and at spatial infinity r* = r*: 

d T ip D>J = DjI , 
d T 4> D ,i = ~d T ip Da , 

5 t^d,i = { 9 r*4>n)j > 



(69) 
(70) 
(71) 



(72) 
(73) 
(74) 



where the subindex 1 refers to the first subdomain and 
D to the last one. 



B. Evolution of the characteristic fields 

An alternative to the penalty method is to stick to 
hyperbolic methods and use the symmetric hyperbolic 
structure of our system in order to communicate the 
different subdomains by passing the characteristic fields 
through the corresponding characteristics. This process 
can be illustrated with the use of Figure [2] Let us con- 
sider the interface between two arbitrary subdomains, say 
Q A and f2 A +i- The boundary nodes are r* a = rj B and 
r* +1 L = . Then, looking at the representation of 
the propagation of the characteristic fields, U tm and V e,m 
in Figure [2| we can deduce that at r* — r* R we only need 
to evolve TEqs. (49 1 for (ip A N , U AN ), whereas V A N is just 



obtained by copying the value of V A+1 from the other 



subdomain 
evolve Eqs 



In contrast, at r* = r 



7 A+l,0> ^A+l,o)l 



A+1,0 
A+1,L 



we only need to 



whereas Y, 



A+1,0 



is 



(491 for (ip 

obtained by copying the value of U A N from the other sub- 
domain. In the case where the particle is at the interface 
of these two subdomains (and hence they are dynamical 
subdomains), i.e. A — P, we cannot just copy the values 
of the characteristic fields, but instead we have to enforce 
the junction conditions of Eqs. (29 1 and (30 1. That is, 
we obtain V p N by doing: 



Vr 



p+1,0 



V 



ip,p+i 



and U p+10 by 



"p+1,0 — ^P,N 



'U] 



(75) 



(76) 



The practical implementation of this technique re- 
quires to treat the different values of the fields at the 
different collocation points in a particular order. In what 
follows we outline an algorithm to do so and the associ- 
ated equations. Step by step, the procedure is the fol- 
lowing: (i) For each subdomain f2 A (A = 1, . . . ,D), we 
evolve the variable U m at the inner points and at the 
right boundary node, r* = r* N , according to the fol- 
lowing evolution equation (again we drop the harmonic 
indices £ and m): 



1 



ip t x A \ 



(d r ,X A \ 



(77) 



where (d r ,X A ) i and (d t X A ) i are given in Eqs. (45 1 
and (46 1 respectively, and i = 1, . . . , N. The next step is: 
(ii) For each subdomain f2 A (A = 1, . . . , D), we evolve the 
variable V 171 at the inner points and at the left bound- 
ary node, r* = r* , according to the following evolution 
equation: 



d T V AA = 



1 - 



(d t X A 



(d r ,V A 



vri>* 



(78) 



where i = 0, . . . , N — 1. (iii) We evolve the variable V im 
at all the right boundary nodes according to Eq. ( 75 1 . For 



the domains without the particle we just simply pass the 
value to the right. Here, we need to make an important 



10 



remark about the practical implementation (in a numer- 
ical code) of this step. For the evolution of the time- 
dependent variables of our system of ODEs, N A i (T), 
we use two different methods. One implements directly 
Eq. ( 75 1 . The other one, uses the method of lines (see, 



e.g. |70|). so that it incorporates the boundary matching 
conditions into the evaluation of the right-hand side of 
the evolution equations. Then, to insert this step into 
the method of lines of our numerical algorithm we need 
to use the derivative of Eq. ( 75 1 rather than this equation 
itself. That is, we need to use: 



dV % 



p+i,o 



d[vr P - 



dT 



dT 



dT 



(79) 



This means that we also need the expressions for the 
derivatives of the jumps in the variables U and V . 
They can be directly derived from the expressions of 
the jumps [Eqs. ( |29"| and (30l] and the geodesic equa- 
tions (see Appendix|A[) . We give the resulting expressions 
in Appendix [c] in particular d[V} F p P+1 /dT is given m 
Eq. flC^ . 

(iv) We evolve the variable U lm at all the left boundary 
nodes according to Eq. (76 1. Again, for domains that 



do not have the particle, we just pass the value of this 
variable, to the left in this case. In order to use the 
method of lines, following the previous discussion, we 
need to use the derivative of Eq. (76 1, that is, 



dU„ 



d[U)\ 



dT 



dT 



dT 



(80) 



ip,p+i 



/dT is given in Eq. (CI I 



where d \ U 

(v) We evolve the characteristic variables U tm and V im 
at the global boundaries. At the horizon, 7-* = r* , based 
on the method of lines, we use the equations: 



0. 



d T v Q>u = (d r ,vg H -v,°' H v>o,H- 

Then, at spatial infinity r* = r*: 

-(9 r »£U-V, D 'V D , r , 



d T U D1 -- 

a T v Dil = o. 



(81) 
(82) 

(83) 
(84) 



(vi) The last step consists in evolving the variable ip im 
at the collocation points of all domains according to the 
equation 



(u . - v ■) 

2 (0 r ,X A ) i 1 A ' J A ' t} ' 

(85) 



where i = 0, . . . , N, and A = 1, . . . , D. 

V. RESULTS FROM THE SIMULATIONS 

In this section, we present the main results from the 
actual numerical calculations. These calculations have 



been done using implementations of the PwP formulation 
in combination with the multidomain scheme described 
above. The details of the numerical codes that we have 
written for the simulations are essentially the same as in 
the case of Paper I, so we refer the reader to this paper 
for the computational details. Here, we focus on the new 
aspects of the implementation. 

The results that we present refer to three different 
types of eccentric orbits (see Paper I for the case of 
circular orbits). These three types of orbits, in terms 
of the eccentricity and semilatus rectum, are: (i) (e,p) 
= (0.1,6.3); (ii) (e,p) = (0.3,6.7); and (iii) (e,p) = 
(0.5, 7.1). We show pieces of these orbits in Figure [3] 

Before introducing the results, it is worth making a 
comment on the numerical implementation of the meth- 
ods and techniques that we described above. In contrast 
with the work done for Paper I, where we only used the 
penalty method to communicate the subdomains, here we 
have implemented in our numerical codes, two different 
techniques: The penalty method, described in Sec. |IV A] 
and the direct communication of the characteristic fields, 
described in Sec. |IV 5] In practice, the numerical imple- 
mentation of the two methods is equivalent to having two 
different numerical codes for the evolution of the system. 
For the results presented in the paper we have mainly 
used the penalty formulation. But we have also used 
the characteristic one, the one based on direct commu- 
nication of the characteristic fields, both to check the 
results, and for the numerical evaluation of the self-force 
below. A last comment on the implementation refers 
to the characteristic formulation. For the time evolu- 
tion, like in Paper I, we use a Runge-Kutta 4 (RK4) 
algorithm (see, e.g. [711 172]). Then, as we already dis- 
cussed in Sec |IVB| there are two ways of communicating 
the fields: Either we use the expressions containing the 
jumps in (U em , V ), Eqs. (76 I and (751, which requires 



some engineering in the RK4 algorithm, or we use the 
expressions contai ning the deri vatives of the jumps in 
(U lm ,V em ), Eqs. (76) and {75}, which can be adapted 
naturally with the method of lines. This second op- 
tion requires to impose initially the values of the jumps, 
Eqs. (CI I and (C2|, since during the evolution the only 



input about them is the information on their derivatives. 
Therefore, in this case, we cannot just use the zero initial 
data that we used in Paper I, and that can be used for 
the other cases presented here, namely U — U(t = t ,r) 
= {tp em (t ,r),(f) im (t ,r),ip lm (t ,r)) = 0, but a modifi- 
cation of it, at least at the nodes where the particle is 
located. All these options have been implemented in our 
numerical codes. 

The first results that we present refer to the valida- 
tion of the code. In Paper I we already presented several 
tests. Here we show a multiple convergence plot (for the 
different orbits) that shows that our PwP scheme indeed 
provide solutions smooth enough to preserve the a priori 
exponential convergence of the PSC method. In Figure [4] 
we have plots for the evolution of the harmonic mode 
(£,m) = (2,2), and the different orbits considered in 
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Orbits around the MBH 





-15 -10 
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10 15 



FIG. 3: Eccentric orbits considered in the numerical evolutions presented in this paper. From left to right we have: (i) (e,p) = 
(0.1,6.3); (ii) (e,p) = (0.3,6.7); and (iii) (e,p) = (0.5,7.1). These orbits have been integrated using the equations of Appendix[A) 
in particular Eqs. (A5l and (A6l. 



Figure [3j showing an estimation of the truncation error, 
[where the a N is the last spectral coefficient of 



10 l"jVI 



log 

the spectral approximation of Eq. (38l], with respect to 
the number of collocation points, N. For these plots we 
use the information coming from the subdomain to the 
right of the particle, i.e. fl r +i, so that the particle is 
located at the left boundary node of this subdomain. As 
one can see, the truncation error, estimated as the abso- 
lute value of the last spectral coefficient, \a N \, decreases 
exponentially with the number of collocation points, as 
expected in the PSC method for smooth solutions. 

The key ingredient of our PwP scheme is to put the 
particle at the interface between two subdomains, so that 
it does not appear explicitly as a distributional term in 
the field equations, but as a boundary condition between 
the solution at adjoining subdomains. This provides 
smooth numerical solutions that satisfy the junction con- 



ditions dictated by the jumps [see Eqs. (23 1 and (24 1 or 



Eqs. (29 1 and ( 30 1 ] . We have performed simulations that 
use zero initial data (or initial data that is zero every- 
where except at the particle nodes), that is 

i> im (t ,r*) = ^ m (t ,r*) = <p tm (t ,r*) = 0. (86) 

As a consequence of this, the numerical evolution pro- 
duces an initial unphysical burst. We have to wait until 
this unphysical feature propagates away in order to an- 
alyze the solution and obtain physically relevant results. 
In Figure [5] we show snapshots of the evolution of the 
(I, m) = (2, 2) harmonic mode, including details near the 
particle location, that illustrate the ability of our method 
to capture the structure of the solution near the parti- 
cle, in particular the ability of resolving the jumps in the 
time and radial derivative of the field variables. 

The main goal of this paper is the computation of the 
self-force, hence the last two plots that we present show 
results on the estimation of the regular field. From the 



values of the components of the gradient of the regular 
field, V Q $ R , the self- force can be directly computed by 
using Eq. ([6]). We have carried out evolutions of the 
orbits shown in Figure [3] and we have computed the 
values of the components of the gradient of the regu- 
lar field along the evolution. This means that we have 
evolved all the necessary modes of the retarded field up 
to a maximum £, say £ majc . To compute the total num- 
ber of evolutions needed we have to take into account 
the following facts: (i) We need to compute the real 
and imaginary parts [the source term in the field equa- 
tions ([8| is complex as shown in Eq. ( p~0~[ )] for each (£, m) 
harmonic component of 3r m ; (ii) the real character of 
the scalar field implies that for each (£, m) the relation 
<jj£-m _ ^_^m^hn ^ we ^ Q n0 ^ neec j £ compute the 

modes with m < 0) holds; and (iii) For any {£, m), we 
have 



Y tm (ir/2,(p) = for l + m odd, 



(87) 



and hence, the jumps for these modes are exactly zero 



[see Eq. ([18) and Eqs. (|2^,([24| and Q,([30|] and we do 

not need to compute them. Therefore, the total number 
of evolutions we need to make in terms of £ m „^ is 



N 



(Cax + 1) (Ca* + 2) 



Once we have evolved the equations for the retarded field 
we have to regularize it, I by £, by computing the compo- 
nents of the gradient of the singular field at the particle 
location, either using Eq. ( 13 1 or Eq. (B16l, and sub- 
tracting it from the gradient of the retarded field. In 
Figure [6] we show the results of the computations of the 
components of the gradient of the regularized field, which 
have been obtained by choosing £ max — 17, for which we 
needed to perform 171 evolutions of the retarded field 
equations. The particular multidomain framework that 
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est radial location to the pericenter radius [see Eq. (A2|] 
available in our evolutions. In Table H] we show these 
values for the three orbits of Figure [3] In order to have 
some reference to compare these values, we also provide 
the values for the last stable circular orbit (r = 6 M) 
obtained, using frequency-domain methods, in |61J . Our 
numerical values have been obtained by using a multido- 
main framework consisting of D — 80 subdomains and 
N = 50 collocation points. Here, we need to make an 
important comment: In Paper I, where the evolutions 
were always done using the penalty method, we always 
quoted two values for each component of the gradient of 
the regularized field. The reason is that in our framework 
we have two different nodes for the particle (one that be- 
longs to the subdomain to the left of it, fl P , and one 
that belongs to the subdomain to the right of it, il P+ i). 
Since the equations are asymmetric with respect to the 
particle location, the evolution via penalty method only 
provides values of the variables at the particle location 
that are approximately the same but not identical. The 
advantage of this is that it provides a criterium to de- 
cide what is the maximum accuracy of the calculations. 
In contrast, the new method of evolution proposed in 
this paper, the direct communication of the characteris- 
tic fields, provides, by construction, identical values at 
both nodes associated with the particle (excepting in the 
case in which we use spectral filters to eliminate spurious 
high-frequency modes in the solution, which is not what 
we have done for the results of this paper) . Here, we have 
presented values obtained with the formulation that uses 
the direct communication of characteristic fields (which 
provides values that are the same, up to machine preci- 
sion, at the two particle nodes), but we have checked that 
they are consistent with the values obtained by using the 
penalty method. 



VI. CONCLUSIONS AND DISCUSSION 



FIG. 4: We show the dependence of the truncation error, 
as estimated from the quantity log 10 \ a N \, on the number of 
collocation points, TV, for evolutions of the harmonic mode £ = 
m = 2 of the field variable tp . From top to bottom, we show 
the truncation error for: (i) (e,p) — (0.1,6.3); (ii) (e,p) = 
(0.3,6.7); and (iii) (e,p) = (0.5, 7.1). The plot corresponds to 
the solution on the subdomain to the right of the particle. The 
tortoise radial coordinate size of this subdomain, i.e. |7"p +1N — 
r*(i)\, is in the range 20 - 40 M. The good fit of the data to 
a straight line confirms the exponential convergence of the 
numerical computations. 



we employed for these computations consists of D = 10 
subdomains and N = 100 collocation points per subdo- 
main. 

Finally, we also want to show some numerical values of 
the components of the gradient of the regularized field. 
To that end, we have chosen to evaluate them at the near- 



In this paper we have extended the time-domain tech- 
niques introduced in Paper I for the computation of the 
self-force from the case of circular orbits to generic or- 
bits, what we have called the PwP formulation. The key 
ingredient of this formulation is to consider a multido- 
main framework in which the particle is always located 
at the interface of two subdomains, even in the case of 
eccentric orbits. In this way, the equations that we have 
solved at each subdomain are homogeneous wave-type 
equations for the retarded field, without distributional 
sources. Then, the solutions of these equations are suf- 
ficiently smooth to provide a good convergence of the 
numerical method involved, in our case the PSC method 
(which provides exponential convergence for smooth so- 
lutions). The way in which we keep the particle at a fixed 
node is by using a time-dependent mapping between the 
physical domain (the one described by the Schwarzschild 
radial coordinate r) and the spectral domain (the one in 
which the calculations of the PSC method take place). 
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r* [M] r * [M] r* [M] 

FIG. 5: Snapshots of the evolution of a scalar charged particle in eccentric orbits around a non-rotating MBH. The orbital 
parameters are (e,p) = (0.5,7.1). The simulations used D = 10 subdomains and N = 100 collocation points per subdomain. 
They show the evolution of the variables ip lm (left), cj> em = d t ip em (center), and cp em = d r -^ tm (right), for the mode I = m = 2. 
In particular, they show (see the boxes inside the plots) how the jumps in the time and radial derivatives of the scalar field 
$' m are resolved in the multidomain computational framework that implements the PwP scheme. 




Time [M] Time [M] 



FIG. 6: Evolution of the components of the gradient of the regularized field, V Q $ R , for a scalar charged particle in eccentric 
orbits around a non-rotating MBH. From left to right, the orbital parameters of the orbits are: (i) (e,p) = (0.1,6.3); (ii) (e,p) 
= (0.3, 6.7); and (iii) (e,p) = (0.5, 7.1). For each orbit (frame), the solid line represents the evolution of the dimensionless time 
component, =y- <3>f ; the dashed line represents the evolution of the dimensionless radial component, $ R ; and the dot-dashed 

line represents the evolution of the dimensionless azimuthal component, ^j- 4> R . The numerical setup for these calculations used 
D = 10 subdomains and N = 100 collocation points per subdomain. 



TABLE I: Numerical values of the components of the gradient of the regularized field at the pericenter radius. We present our 
estimations for the orbits shown in Figure 3] (i) (e,p) = (0.1,6.3) [Column 3]; (ii) (e,p) — (0.3,6.7) [Column 4]; and (iii) (e,p) 
= (0.5, 7.1) [Column 5[. For reference, we have included the values for circular orbits at the last stable circular orbit (r = 6M) 
obtained, using frequency-domain methods, in |61| [Column 2[. 





(e,p) = (0.0,6.0) 


(e,p) = (0.1,6.3) 


(e,p) = (0.3,6.7) 


(e,p) = (0.5,7.1) 


^'peri (^0 


6.0 


5.7272727 


5.1538461 


4.7333333 


r p (M) b 


6.0 


5.7272925 


5.1538801 


4.7333989 


q t 


3.609072 ■ 10" 4 


4.5171 ■ 10" 4 


7.6980 ■ 10" 4 


1.2330- 10" 3 


9 


1.67728- 10" 4 


2.1250 • 10" 4 


3.6339 ■ 10" 4 


5.6122 ■ 10~ 4 


q 


-5.304231 ■ 10~ 3 


-6.2040 ■ 10" 3 


-9.0402 • 10" 3 


-1.2685 ■ 10~ 2 



"This is the exact value of the pericenter radius for each orbit. 

'This is the value of the radial coordinate at which the gradient of the regularized field has been evaluated. It is also the nearest value 
of the radial coordinate to the pericenter value available in our numerical evolutions. 
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The details of this mapping are given in Sec. |IV| The 
technical difficulties are then transfered to the commu- 
nication between the different subdomains, in particu- 
lar between the subdomains that surround the particle. 
In this sense, we have extended the use of the penalty 
method to the eccentric case and proposed a new method 
for the communication between subdomains based on 
the transfer of the characteristic fields of the hyperbolic 
system of equations. These two methods are different 
enough that their numerical implementation can be con- 
sidered to lead to different numerical codes, which is good 
for checking the results of the computations. 

The main advantages of the PwP formulation in com- 
bination with the PSC numerical method were already 
described in Paper I and are still valid for what we have 
presented in this paper, so we do not discuss them fur- 
ther here. Instead, we discuss the possible improvements 
of the methods presented. In Paper I, we already men- 
tioned a number of such improvements: (i) Compactifi- 
cation of the computational domain, in other to improve 
the accuracy of the (global) outgoing boundary condi- 
tions which will allow longer evolutions maintaining the 
accuracy, (ii) Reduction of the time step, by abandon the 
linear mapping between the physical and spectral repre- 
sentations, using a mapping of the type introduced by 
Kosloff and Tal-Ezer [73] • (iii) Richardson extrapolation, 
to improve the estimations of the values of the self-force 
by using our analytical knowledge of the expansions in 
inverse powers of the harmonic number £, like in |47j . 
These improvements, although they have not been im- 
plemented in the present work, can be perfectly applied 
to our framework and have significant potential to im- 
prove the efficiency of the computations. Another im- 
provement that can be even more important than the 
ones already mentioned has to do with the distribution of 
subdomains. As it was discussed in Paper I, it is very im- 
portant to choose appropriately the (tortoise coordinate) 
size of the subdomains in order to resolve accurately the 
different harmonic modes. In this sense, the modes with 
high harmonic number m are more difficult due to the 



dependence of the source term in the original 



e imip p (t) 

equations [see Eqs. Q and ( 1 1] . This means that the 
higher m is, the higher the frequency of the scalar field 
modes generated is, and hence the shorter the wavelength 
gets. Therefore, these high m modes require more spa- 
tial resolution (also in time). Both in the calculations of 
Paper I and in this paper we have used a unique sub- 
domain structure for all the modes. Since we must set 
the size of the subdomains in such a way that the high 
m modes are properly resolved, this means that we are 
using excessive resources for the low m modes. Thus, the 
optimal strategy would be to adapt the subdomain struc- 
ture to the type of mode we are dealing with. One can 
improve in this way the efficiency and reduce significantly 
the computational resources involved. 

Another way of improving the performance of the nu- 
merical computations is to provide better initial data, or 
at least to reduce the negative effects of the initial burst 



produced when we prescribe zero initial data [Eq. ( 86 1] . 
This is important since as it has been suggested in |74|. 
certain components of this unphysical burst can persist 
in time, which can deteriorate significantly the precision 
of the calculations. A possible way of alleviating the sit- 
uation is to switch on the particle during a certain time 
scale, by multiplying the source terms (in our case the 
jumps) by a convenient function of time (starting at zero 
initially and approaching unity) . We have been exploring 
this possibility, that turns out to be critical in the imple- 
mentation that makes use of Eqs. (73) and (74), and we 
will report on this elsewhere. 

To finish, the main goal of the formulation presented in 
this paper is to develop an accurate and efficient method 
to compute the self-force in situations of physical inter- 
est. In particular, for systems of interest for the future 
observatory LISA. This means to use these techniques in 
the gravitational case and for spinning MBHs. In this 
sense, we have to mention that while it is not difficult to 
transfer the techniques discussed here to the gravitational 
case, to do the same with the case of a spinning black hole 
may require new technical improvements which we will 
the subject of future investigations. 
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Appendix A: Particle motion 

In this section we summarize the form of the geodesic 
equations for orbits around a Schwarzschild MBH that 
we use in this work. Using the constants of motion and 
the normalization condition of the particle's velocity u^, 
g flv u ,J 'u v = — 1, we can reduce the geodesic equations to 
first-order ODEs for (t p ,r p ,ip p ). To avoid the turning 
points in the integration of the resulting equation for the 
radial coordinate, we use an angle variable that is related 
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to the radial coordinate via 



pM 



1 



ecos x p 



(Al) 



where p is the dimensionless semilatus rectum and e is the 
eccentricity. From this expression, the minimum (peri- 
center, r peri ) and maximum (apocenter, r apo ) values of 
the radial coordinate, the turning points, are 

pM pM 



1 



1 



(A2) 



Since the spacetime is static and spherically symmetric, 
the geodesies have the following constants of motion: En- 



ergy E p = —u t and angular momentum L p = u$, which 
can be written in terms of the orbital parameters (p, e) 
as 



Ei = 



T 2 



(p-2-2e)(p-2 + 2e) 



p(p 
p 2 M 2 
p — 3 — e ; 



3-e 2 ) 



The equations for (x p (t),ip p (t)) are 



(A3) 
(A4) 



dx P 
dt 

defp 
dt 



(p - 2 - 2e cos x p ) y/p-6- 2ecosx p (l + e cos x p ) 2 
MpV(p-2) 2 -4e 2 

(p - 2 - 2e cos x P ) (1 + e cos x P ) 2 
Mp 3 / V(P - 2 ) 2 " 4e2 



(A5) 
(A6) 



Two additional expressions of interest for this paper 
are the first and second time derivatives of the tor- 
toise coordinate, f* and r*. They can be obtained from 
the geodesic equations doing appropriate manipulations. 
The result is 




and 



[p 
r p E 2 



M 



T 2 



3M 




(A7) 



(A8) 



Appendix B: Structure of the Singular Field 

In the description of the problem in Sec. [II] it was made 
clear that the computation of the self-force requires the 
regularization of the retarded scalar field. This can be 
done in the framework of the mode-sum scheme, where 
we need to compute both the harmonic components of the 
retarded and singular parts of the scalar field. The singu- 
lar field near the particle has a particular dependence on 
the harmonic number £, which is fully determined by the 
specification of a set of regularization parameters (which 
can be derived by analyzing the field equations in a world- 
tube close to the particle's worldline). In this Appendix 
we summarize the expressions for these parameters that 
we have used in this work (see [461 14"51 E3J E3])- 

We start with the prescription introduced by Haas and 
Poisson j3S], which is the only one that includes the reg- 
ularization parameter that goes with the l~ 2 (for large i) 



behaviour of the singular scalar field, namely D a . This 
term, although it is not necessary for the computation 
of the regularized field, it is useful in order to accel- 
erate the convergence of it as we increase the value of 
l max (the maximum £ harmonic number that contributes 
to the computation of the regularized field, and hence 
to the self-force). Following [55] we decompose V„$ 
in terms of a complex tetrad e?^(ai), with the index 

(m) = {(0), (+), ( — ), (3)} labeling the different compo- 
nents of this tetrad 



'(0) 



'(3) 



; (±) 



1 

77 



0,0,0 



0,^0080,-^,0 



0, y/fwaOe^ 



(Bl) 
(B2) 
,(B3) 



The components of the gradient of the retarded field in 



this tetrad, $ 



$ e? , 



are scalars. Thus, we can 



expand $/ in scalar spherical harmonics: 



$ (m)(^) =EE *W){t,r)Y tm {9,<p) . (B4) 



£=0 m= 



Then, the components of the gradient of the regular field 



[Eq. ( 11 1] are given by 



$ s = 



" $ (M) 



(B5) 



where, $^ and ^f^, are the tetrad projections of the 
gradient of the retarded and singular fields. Then, the 
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relation between the tetrad and coordinate components 
of the gradient of the regular field can be written in the 
form: 



The components A,^ of the regularization parameters 
are given by 



®(0) 



*(+) _ 



1 

77 

-a 



*(-) ~ 
*(3) - 



(-) 

) + /3| m (V+ lm ) 



& tm W-lm 



(B6) 

,(B7) 
(B8) 
(B9) 



where dfe, /3fc, orfe, /fe, a^ m , and /3| m are linear 
operators given by the following expressions 



*(+) 



£-1 



' (1 + to- 1)(1 + to) f r-d_ 

(2£-l)(2£+l) y } dr r 



(BIO) 



^(0) = 



.4 



e E r 
p p p 



(+) 



^(3) = 



A, 



(B17) 
(B18) 



= -e Ev ^ V 
V \Hp ( t p ^p) 
0, (B19) 
where we recall that a dot means differentiation with re- 
spect to the coordinate time t and e p is 



e p = sign(r - r p ) 



1 if r > r p 
— 1 if r < r„ 



(B20) 



lm (i-m-l)(i-m) 
<"> V (2£-l)(2£ + l) 



/3^m 



03 = 



4-^1 (BID 

or r 



.,., = „ fc^gj^ r^-^'i. (Bi2) 



(2f- 1)(2£ + 1) V or 





— TO - 


f 1)(*- 


- 771 - 


-2) 




(2£- 


f l)(2f 


+ 3) 






+ TO - 


f 1)(^4 


- TO - 


-2) 




(2£- 


f 1)(2* 


+ 3) 




/"- 


- TO -f- 


w+ 


TO + 






(2^4 




-3) 





>r r 



B13) 



'4 + ^W) 



9r r y 



or r 



(915) 



On the other hand, the tetrad projections of the sin- 
gular field can be written [in a similar way as in the 
coordinate case of Eq. ( 13 1] as 



(a0 



(B16) 



The components B,^ of the regularization parameters 



are 



E 2 r f 
p p p 



2[fp(r 2 p + Ll)] 



3/2 V V 



^i/2-2F_ 1/2 , (B21) 



where and F_ ± ^ 2 are objects defined in terms of 

hypergeometric functions |75| as 



F -^- F{ 2'2 ;1; r2+ P L2)' 



(B22) 



(B23) 



and 



B (+) = B ( _) = e** - »[«(+)]) , (B24) 



(+)J 



'(+)J 



where 3?(z) and 5(z) denote the real and imaginary parts 
of the complex number z. In the case of , these parts 
are given by 



(+)J 



E 2 p 



2fp /2 (r 2 P + Ll) 



3/2 



(2F_ X 



/2 - F l/2 



^-1/2-2(^-1)^ 



/2 



[*<+)] = 



(2-v%) E p r p 



2L p / p 3/2 v /(r2+L2) 



^1/2 ~ ^-1 



/2 



(B25) 



(B26) 
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Moreover, we have Finally, the components of the regularization pa- 

rameters are 

B {3) = 0. (B27) 



and 



C M =0. (B28) 



+ Mr 2 p (l6 F 1/2 - 30 F_ 1/2 ) - r p L 2 p (3 F 1/2 - 7 F_ 1/2 )+ML 2 p (42 F 1/2 - 114F_ 1/2 ) 
M£4 , s 36ML? 



+ 



- (18 F 1/2 - 104 F_ 1/2 J - —^F_ 1/2 }> , (B29) 



%) = A-) = e ^ p - 1 . ( B3 °) with 



Elfl if, , , o 26ML1 9ML«\ 
+ ^77 \ 4 3r P + 6Mr P - L P r P + 31Mi P + + T 2 - F -i/2 

- \7r 3 p + YIMvl - L 2 p r p + A&MI 2 + \ F 1/2 j 

JL if Ll 2QML1 22MLt 8ML 6 V \ 

16(r2 + L2)3/2 ^ P Tp r 2 r 4 r 6 J I 

f QLl \2MLl 4MLt\ ) 

1 "^-^W^l (B3D 



/3 P P 



£Lr_ (/ „ 2 10ML 2 , 



l6L p f p (rj + Liy 



p 



10MLI 29MLt UML 6 V \ 



- 2 ^ + Mr p + 5L; + + + j F_ 1/2 

E v r v i f , „ 5M£2 2MLl\ 

+ -^TT^ 2 K"^-ir + ir T-V2 



2r^ - 2Mr p + 5L 2 p - F 1/2 , (B32) 
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and 



r> (3) = o. 



(B33) 



On the other hand, the regularization parameters as- 
sociated with the coordinate expression of the singular 
field [see Eq. ( 13 1] , given in |46l IBUl IM] . are given by 



-4, 



E r 
p p 



(B34) 



TV IF 

f p 



C = . 



F 



-1/2 



(B39) 



(B40) 



E„ 



'fp(r 2 p + L 2 p)' 



In this case, the regularization parameters D a have not 
^ ' been computed. 



^ = 0, 



(B36) 



E 



B, 



' 2 r r ( F-, 
p p p \ l 



/2 2 ^-l/2 



2f p (r 2 +L 2 p ) 



3/2 



(B37) Appendix C: Time Derivatives of the Jumps in the 

Characteristic Variables 



E 2 r 
p p 



(r 2 p -2f 2 )F 1/2 + (r 2 -f 2 )F_ 1 



/2 



3/2 



The time derivatives of the jumps of the characteristic 
,(B38) fields, are given by the following expressions: 



d[U] p 
dT 



Airq f. p 

-^p^p 1 ^p 

47rg f p 

E r 2 1 - r* 
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E p r P 
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47rg f p 
E p r 2 p 1 + r; 
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E r 
^p' p 

m Vp m Y ^] 

EpVp 



(C2) 



where 3t[Y ■ ] and SJ[Y ,m ] are the real and imaginary r* are given by Eqs. (A7l and (A8l respectively, 
parts of the scalar spherical harmonic Y em , and r* and 
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